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Abstract 

We solve exactly the problem of a finite slab receiving an isotropic radiation on 
one side and no radiation on the other side. This problem - to be more precise 
the calculation of the source function within the slab - was first formulated by K. 
Schwarzschild in 1914. We first solve it for unspecified albedos and optical thick- 
nesses of the atmosphere, in particular for an albedo very close to 1 and a very large 
optical thickness in view of some astrophysical applications. Then we focus on the 
conservative case (albedo = 1), which is of great interest for the modeling of grey 
atmospheres in radiative equilibrium. Ten-figure tables of the conservative source 
function are given. From the analytical expression of this function, we deduce 1) 
a simple relation between the effective temperature of a grey atmosphere in radia- 
tive equilibrium and the temperature of the black body that irradiates it, 2) the 
temperature at any point of the atmosphere when it is in local thermodynamical 
equilibrium. This temperature distribution is the counterpart, for a finite slab, of 
Hopf's distribution in a half-space. Its graphical representation is given for various 
optical thicknesses of the atmosphere. 

Key words: Radiative transfer; Plane-parallel geometry; Isotropic scattering; Grey 
atmosphere; Radiative equilibrium; Local thermodynamical equilibrium; 
Temperature distribution. 
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1 Introduction 



In a paper considered as foundational in transfer theory, K. Schwarzschild [1] 
modeled the solar atmosphere as a finite slab irradiated on one side by black 
body radiation. He set up, for the first time, the complete radiation transfer 
equation - including the scattering integral - and inferred the following integral 
equation for the source function S: 

S(r)=S (r) + y^E 1 (\r-r'\)S(r')dT', (1) 

with 

S (t) = (1 - a)B*(r) + ^B(T b )E 2 (b - r). (2) 

The first term on the right-hand side of Eq. (1) describes the contribution 
of the (internal and external) sources, and the integral term is the scattering 
term, a G [0, 1] is the albedo (supposed constant) for single scattering by a 
volume element, b > the optical thickness of the atmosphere, and r G [0, b] 
is the optical depth variable. Eq. (2) shows that the source term S includes 
the thermal emission B* of the atmosphere, as well as the contribution of 
the external sources. In our model, they will be prescribed as black body 
radiation of temperature T b , incident on the surface r = b. B(T b ) is the Planck 
function at the temperature T b . There is no radiation incident on the other 
boundary surface r = 0. The scattering of light is supposed monochromatic 
and isotropic, which assumption gives rise to the exponential integral functions 

EJt) = f 1 exp(-r/u)u n - 2 du (n > 1) (3) 
Jo 

in the kernel and the free term of Eq. (1). Frequency dependence of quantities 
will not be made explicit. 

The integral form of the radiation transfer equation was established before the 
publication of Schwarzschild's paper, by Lommel [2], Chwolson [3] and King 
[4]. In these articles, Eq. (1) with S (t) = exp(— r/fi) (/i > 0) was derived 
directly, with no reference to the transfer equation. Integral equations of the 
form (1) have been studied intensively since these pioneering works. Three 
general methods for solving them, which apply to the particular source term 
(2), are described in [5]. Since the problem (l)-(2) is linear, its solution can 
be written as 

S(r) = (1 - a) f R(a, b, r, r')B\r') dr> + B(T b )Z (a, 6, 6 - r), (4) 
Jo 

where R denotes the resolvent kernel of the integral equation (1) and £ the 
solution to the following integral equation: 

Ua, b, r) = a -E 2 {r) + | £ E^r - r'|)£ (a, b, r') dr'. (5) 
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In view of evaluating the integral on the right-hand side of Eq. (4), the 
model under investigation must be further specified to make explicit the r- 
dependence of the function B*. When the distribution of internal sources is 
uniform in the slab, B* is independent of r and then the thermal emission of 
the atmosphere is characterized by the function 



In the conservative case (a = 1), both functions £ an d Q were studied in detail 
by Heaslet and Warming [6]- [7] in view of some applications in heat transport 
theory Their functions O(t) and B s (r) coincide with our functions £o(l,b, r) 
and (1/4)Q(1, b, r), respectively. In the non-conservative case (0 < a < 1), 
the function £ is briefly evoked in Section 8 of [8], while the literature on the 
function Q is mentioned in a recent article devoted to this function [9]. 
In what follows, we will be interested exclusively in the second term on the 
right-hand side of Eq. (4), which describes the contribution of the external 
sources. The exact expression for the £ -function is given for an unspecified 
albedo a, including a — 1. This particular case is of great importance for 
our understanding of grey atmospheres in radiative equilibrium. Actually, it 
is known [10] that the solution to Schwarzschild's problem with a = 1 is the 
source function, integrated over frequency, of a grey atmosphere in radiative 
equilibrium. It thus yields the temperature of the atmosphere if the latter 
may be assumed to be in local thermodynamical equilibrium. This particular 
problem has been solved exactly in a half-space, that is, for a = 1, b = +oo 
and Sq = (since the sources are at infinity). The solution is due to Hopf [10]. 
In a finite slab (b < +oo), the conservative problem (l)-(2) is equivalent to (5) 
with a — 1. Yamamoto [11] failed to solve it exactly, as remarked by Sobolev 
[12]. Heaslet and Warming's solution [6]-[7] is approximate, in contrast with 
that given by Kriese and Siewert [13], who used the Case method. Four- figure 
tables of the conservative functions £ an d Q are given in [7] , in good agree- 
ment with the six- figure tables of [13]. 

The outline of this article is as follows: in Section 2, a recent work dealing with 
a problem connected with Schwarzschild's problem is summarized. Sections 3 
to 5 are devoted to the calculation of the £ -function in the general case, in a 
conservative atmosphere and in a semi-infinite atmosphere, respectively. The 
numerical evaluation of the £o~function is tackled in Section 6, which contains 
a table of this function for a = 1 and different values of b. In Section 7, an ap- 
plication of the conservative case to grey atmospheres in radiative equilibrium 
is investigated. From the exact expression of the source function integrated 
over frequency, we deduce a simple relation between the effective temperature 
of the atmosphere and the temperature of the black body that irradiates it. 
We also give the temperature of the atmosphere when it is in local thermo- 
dynamical equilibrium, a calculation which generalizes Hopf 's calculation in a 
half-space to a finite slab. 




b 
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2 The standard transfer problem in a stellar atmosphere: a sum- 
mary 



In the simplest model atmosphere that can be conceived - a homogeneous 
and isothermal slab in local thermodynamical equilibrium with no incident 
radiation - the key function is the function Q, which is the solution to the 
following integral equation 

Q(a, b, t) = 1 + £ f E x {\t - r'\)Q(a, b, r') dr>. (7) 
2 Jo 

It is clear that this function is symmetric about r = b/2, i.e., Q(a,b,r) = 
Q(a, b,b — t). A detailed, recent study of this classical auxiliary function can 
be found in Ref. [9], which we summarize hereafter. The starting point is 
Sobolev's expression [14] 

Q(a, b, t) = ip(a, b, b)[ip(a, b, r) + ip(a, b,b - r) - ip(a, b, &)], (8) 

which involves the function 

^M,t) = 1 + f T 0(a, b, r') dr'. (9) 
J o 

Here, = 0(a, 6, r) is the resolvent function of the problem (1), which solves 
it for the particular free term Sq(t) = (a/2)Ei(r). From the well-known an- 
alytical expression of this function, we derived in [9] the following expression 
of the -^-function: 

a (a,b) 

ip{a,b,T) = — 

l — a 

-R(a){[l-^ x (a,b,l/k(a))]exp[-k(a)r]+^ Y (a,b,l/k(a))exp[-k(a)(b-r)]} 

~o / 9(a,u) {[1 - £ x (a,b,u)] exp(-r/u) + £y(a, 6, w) exp[-(6 - t)/u]} du. 
2 Jo 

(10) 

Let us briefly point out the definition of the coefficients and functions appear- 
ing in this expression, referring to [9] for details: 

1) For < a < 1, the coefficient k(a) is the unique root, in ]0, 1[, of the 
characteristic equation 



2k{a) 



1 + k(a) 
1 - k(a) 



0, (11) 



so that k(0) = 1 and k(l) = by continuity. 
2) We have 



1 — k (a) . . . 
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3) For any x ^ —I 

£x(a,b,x) = ^-x f X(a,b,u)-^—, (13) 
2 Jo u + x 

fr(a,b,x) = ^x [ Y(a,b,u)-^-, (14) 
2 Jo u + x 

where X and Y are the classical functions of finite, plane-parallel media. Note 

that the integrals are Cauchy principal values when x e] — 1, 0[. 

4) We introduce the moments of order n > of the X- and F-functions 



,(a,b) = [ X(a,b,u)u n du , f3 n (a, b) — [ Y(a,b,u)u n du, (15) 



and we set 

a (a, b) = 1 - ^a Q (a, b) , j3 (a, b) = |/5 (a, b). (16) 

These coefficients are not independent, since they satisfy 

[a (a, b) - f3 (a, b)} [a (a, b) + (3 (a, b)} = 1 - a. (17) 
5) Finally, the integral in the right-hand side of Eq. (10) contains the function 

9(tt - M) = r>, M ) + V/2)H 2 (°£«<D. («) 

where 

r(a,«) = l--«ln(— ). (19) 
The boundary values of the -^-function are [9] 

^(a,6,0) = l, (20) 
il>(a,b,b) = X(a,b,oo), (21) 



where 



X(a, b, oo) = \a {a, b) - (3 {a, 6)1 = ^ (22) 



1 -a 



a (a,6) +/3 (a, 6) 



is the value at infinity of the X-function; whence the boundary values of the 
Q-function: 

Q(a,6,0) = Q(a,6,6)=A'(a,6,oo). (23) 

The expression (10) for ip(a, b, r) contains the indeterminate form oo — oo when 
a — > 1, since a (a, b) is bounded and R(a) diverges like 1/(1 — a). In [9], the 
functions F± were introduced to circumvent this difficulty Another approach 
would have been to consider that the function ip is the basic function of the 
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problem and to remove the indetermination analytically. This can be done by 
remarking that if) (a, b, 0) = 1, and thus we have 



1 =^3^ " R(a){[l - Zx(a, b, l/k{a))\ + £ Y (a, b, l/k(a)) exp[-A;(a)&]} 

a rl 

-«/ 9(a,u) {[1 -£ x (a,b,u)] + fr(a,b,u) exp(-b/u)} du. (24) 
2 Jo 

Subtracting Eq. (24) from Eq. (10), one obtains the following alternative ex- 
pression for 

ip(a, b, r) = 1 + k(a)R(a)r{[l - ^ x {a, b, l/k{a))Yf[\, k(a)r] 

- £y(a, b, l/fc(a))7*[l, —k{a)r] exp [-k(a)b]} 

+ r a - £ g(a, u) {[1 - fr(a, b, u)] 7 *(l, t/u) 

- fr(a, b, «)y(l, -r/«) exp(-6/«)} — , (25) 

where the function 

7 *(l,x) = -[l-exp(- a ;)] (26) 

x 

belongs to a class of special functions derived from the incomplete gamma 
function: see Ref. [15], Section 6.5. We note that for x > 0, 7 *(1, x) decreases 
from 7*(1, 0) = 1 to 7*(1, +oo) = 0, and that 7*(1, x) ~ l/x — > as x — > +oo. 
As a — > 1, fc(a) — > and 7*(l,±/cr), exp[— fc(a)6] tend to 1. The first term 
in braces on the right-hand side of Eq. (25) tends to 1 — ^x{o-,b, l/k(a)) — 
£y(a, b, l/k(a)), which is proportional to k(a): see Section 4. We thus make 
appear k 2 (a)R(a), which remains bounded as a — > 1. 
It is clear that the function 

4>(a, b, r) = ip(a, b, b - r) - ^(a, b, b) (27) 

can be written in a form similar to (25) by subtracting from both members of 
Eq. (10) their value for t — b. One obtains 

4>(a, b, r) = -k(a)R(a)r {[1 - &(a, b, 1/Jfe(a))]7*[l, -k(a)r] exp[-k(a)b] 

-£ Y (a,b,l/k(a)h*[l,k(a)T]} 

~ T 2J 9 ( a,u ^ V 1 ~ ^x(a,b,u)]-f*(l, -t/u) exp(-6/u) 

-fr(a,M)7 , (l,T/tt)}-. (28) 

J 7/. 



6 



It follows from Eqs. (27) and (20)-(21) that the surface values of the function 

■ip are 

^M,0) = 0, (29) 
4>(a,b,b) = 1 -X(a,b,oo). (30) 

In terms of the functions ip and ip, the expression (8) of the Q-function becomes 

Q(a, b, t) = ip(a, b, b)[ip(a, b, r) + ip(a : b, r)]. (31) 



3 The calculation of £ (a,b,r) 

A natural approach consists in replacing E 2 (t) by its definition E 2 (t) = 
Jq exp(— r/u)du on the right-hand side of Eq. (5). Taking advantage of the 
linearity of this equation, we obtain the solution in the form 

a rl 

£o(a,b,T) = - B(a,b,r,u)du, (32) 

1 Jo 

where the function B is the solution to the following integral equation: 

B(a,b,T,u) = exp(-r/u) + C E 1 (\r - r'\)B(a,b,r' » dr' . (33) 

2 io 

This function is a classical auxiliary function in plane-parallel geometry, which 
coincides with the X- and F-functions on the two boundary surfaces [16] 

B(a,b,0,u) = X(a,b,u) , B(a, b, b, u) = Y(a, b, u). (34) 

Inserting these expressions into Eq.(32), we derive the surface values of the 
£o-function 

£ (a,b,0) = \ {a,b) , £ (a, b, b) = ^f3 (a, b). (35) 

Now derive both members of Eq. (5) with respect to r, remembering that 
-^2( r ) = ~E\(r) ) and note that the solution to Eq. (1) with So (t) = (a/2)£ , 1 (r) 
is the resolvent function 0. Replacing £o( a > b, 0) and £o( a ? b, b) by their expres- 
sions (35), one gets 

e (a,6,r) = ^ /V(|r-r'|)£ (a,&,T')dT' 
I Jo 

- a (a, b^E^r) - p (a, b^E^b - r). (36) 

Hence 

£ (a, b, t) = -a (a, b)<f>(a, b, r) - f3 (a, b)(f)(a, b,b-r), (37) 



7 



and 

£o{a,b,T) = -a (a,b) f <j>(a, b, t 1 ) dr'-/3 (a, b) f (f)(a,b,b-r') dr'+A. (38) 

JO JO 

It follows from Eq. (35) that the constant of integration A is (a/2)a (a : b). 
Using the definitions (9) and (27) of the functions ip and one finally obtains 

£ (a, b, t) = 1 - a (a, b)ip(a, b, r) + f3 (a, 6)V>(a, b, r). (39) 

As a check, the boundary expressions (35) of the £ -function are retrieved by 
means of Eqs. (20)-(21) and (29)-(30). 

Another approach for the calculation of the £ -function consists in substituting 
in Eq. (32) the expression for B(a,b,r,u) given in [17], Eqs. (66)-(68). The 
integration with respect to u can be performed analytically, and the solution 
(39) is derived by transforming the integral along the imaginary axis with the 
help of the residue theorem. 
Introducing the escape probability function [18] 

P{a,b,T) = l-(l-a)Q{a,b,r), (40) 

it follows from Eqs. (39), (27), (21), (22) and (8) that 

£o(a, b, t) + f (a, b,b-r) = P(a, b, r), (41) 

which clarifies the link between the functions Q and £ - I n a conservative 
atmosphere, P(l, 6, r) = 1 and the above equation reduces to the already 
known relation [6] 

£o(1At)+£ (1,M-t) = 1, (42) 

which shows that the conservative £o-function is antisymmetric about the point 
r = 6/2, where it is equal to 1/2. 

The probabilistic interpretation of the relation (42) is clear. We remind the 
reader that (l/47r)S(a, b, r, u)dfl is the probability that a photon at level r 
leaves the slab, after multiple scattering, within an element of solid angle dfl 
around a direction defined by the angle arccosw with respect to the outward 
normal at the top of the atmosphere. From Eq. (32), we deduce that £o(o, b, r) 
is the probability of emergence, through the upper boundary plane, of a photon 
incident on a particle placed at level r. £ (a,b,b — r) is the probability of 
emergence from the lower boundary plane, and the sum of these two functions 
is the probability of emergence P(a, b, r). 



4 The case a — > 1 

The numerical calculation of £,o(a, b, r) with the help of Eqs. (39), (25) and (28) 
is easy as long as the albedo a is not too close to 1 and the optical thickness b 
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is large enough to make sure that the product k(a)b is greater than 30. When 
k(a)b < 30, a difficulty arises when computing the functions ip and ip, since 
the product k(a)R(a), which appears in the right-hand side of Eqs. (25) and 
(28), diverges like l/k(a) as a — > 1. Since it is multiplied by a quantity which 
tends to be proportional to k(a) as a — > 1, an indeterminate form oo x 
appears when a — > 1, k(a) — > 0. To remove it analytically, we write 

a ij - M = i(a + fa)kl> ~ + \{a - #>)(^ + ^) (43) 
with which Eq. (39) becomes 

fo(a, 6, r) = 1 - ^[«o,+(a, 6)^- (a, 6, r ) + «o ,-(a, 6)^+(a, 6, r)], (44) 

where 

a ,±(a,6) = a (a, 6) ± (3 (a,b), (45) 

and 

i/)±(a, b, t) = ^(a, 6, r) ± if) (a, b, r). (46) 
From Eqs. (25), (28) and the identity 

7 *[1, k(a)r] ±7*[1, -fc(a)r] exp[-Jfe(a)6] = 7 *[1, fc(a)r]{l ± exp[-fc(a)(6 - r)]}, 

(47) 

the expression for ip±(a, b, r) is seen to be 

V>±(a, 6, r) = 1 + fc(a)i2(a)r[l - b, l/k{a)) ± £y(a, 6, l/k{a))\ 

x 7 *[l,A;(a)r]{l T exp[-A;(a)(6-r)]} 

+ / ^( a > M )[ 1 -^x(a,6,M)±^y(a,6,M)]7*(l,r/M){lTexp[-(6-r)]} — . 
2 Jo u 

(48) 

The integral term in this expression involves the functions £x an d £y, which 
are difficult to compute on [0, 1]. As in Section 5 of [9], we utilize the classical 
^/-function and the functions (±, first introduced in Ref. [19], which makes it 
possible to calculate the functions X, Y, £x and £y in one step. The algorithm 
to calculate them is summarized in Section 6.3 of [9]. Formulas expressing the 
functions £x and £y are, for x > 0, 

l-£ x (a,b,x) = #£^(C+ + C-)M,aO, (49) 

M«, b, x) = — ^(C+ " C-)(a, 6, (50) 
-H (a, 2 

Consequently, 

1 - £ x (a, b, x) ± £y(a, 6, x) = — - r (±(a, b, x), (51) 
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and the expression (48) for ip±(a, b, r) now reads 

x 7 *[l, fc(a)r]{l =F exp[-k(a)(b - r)]} 

+ ^ / 1 ^(a,«)C ± (a,6,u) 7 *(l,r/«){l=Fexp[-(6-r)]}— }. (52) 
2 jo u J 

Since 1 — exp[— — r)] = k(b — r)7*[l, — r)], we have 

iMa, 6, r) = 1 + r(b - r) { ^^ )) M«)C + («, 6, !/*(<*)) 

x 7 *[l,A;(a)r]7*[l,A;(a)(6-r)] 
+ ^ /^(a !M )C + (a,6, M )7*(l ) r/ M )7ll ! (&-r)H^), (53) 

and 

x 7 *[l, fc(a)r]{l + exp[-A;(a)(& - r)]} 
+ ^ / 1 y(a,«)C-(a,6,«) 7 *(l,r/ U )){l + exp[-(6-r)/ U ]} — ). (54) 

As a — > 1, the coefficients k(a)R(a)/H(a,l/k(a)), k(a)( + (a,b,l/k(a)) and 
C-(o, b, l/k(a)) remain bounded, with the following limits [9] 

k{a)R{a) U = V5, (55) 



H(a,l/k(a)y 



k(a)( + (a,b,l/k(a))\ a=1 = ^A,(l,6), (56) 

C_(a,fe,lA(a))| a =i = ^(«i + /3i)(l,fe). (57) 

Moreover, asymptotic expansions (not given here) of these coefficients for a — > 
1 can be written, which allows one to calculate them for values of the albedo 
very close to 1. It is thus easy to compute the functions tp± with the help of 
Eqs. (53)-(54), even when a — > 1. For a = 1, we have from Eqs. (55)-(57) 

^ + (1,6,t) = 1 + t(6-t){|a,(1,6) 

( / (l,«)C + (l,6,«) 7 J,, (l,r/«)7*[l,(6-r)/«]^|, (58) 
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and 

^_(l,6,r) = l + r{^(a 1 + /3 1 )(l,6) 

+1 / 1 ( ? (l, M )C-(l ) & ) «)7*(l,r/ M )){l + exp[-(6-r)H} — ). (59) 
2 Jo ii J 

On the other hand, the coefficients ao(a,b) and j3 (a,b) remain bounded as 
a — > 1, with limits verifying a (l, 6) + A)(l, b) = 2 [16]. As a result, 

<5 ,-(1,6) = , a 0) +(l>6)=A)(l,6) (60) 
and Eq. (44) reduces for a = 1 to 

e (l,6,r) = 1-^(1, 6)0_(1, 6,r). (61) 
An alternative expression given in Ref. [6] is 

e (l,6,r) = ^[l + A)(l,6)F_(l,6,r)], (62) 

where F_(l, 6, r) = 0(1, 6, 6 — r) — 0(1, 6, r) is the function already introduced 
in Ref. [9]: see Eq. (43), which for a = 1 becomes 

F_(l,6,r) = |(ai + ) 8 1 )(1,6)(6-2T) 

n / 1 -^ 1 4C-(l,6,«){exp(-r/«)-exp[-(6-r)/«]}d«. (63) 
2 Jo H (1, ii) 

To infer Eq. (62) from Eq. (61), just replace ^_(1, 6, r) by ^(1, 6, r) -0(1, 6, r) 
= ^(1, 6, r)— -0(1, b, 6— r)+-0(l, b, b) into Eq. (61), and observe that if)(l, b, r) — 
■0(1,6,6 — r) = — F_(l,6, r) [as seen by comparing the Eqs. (36) and (42) of 
[9]] and 0(1,6,6) = X(l,6,oo) = 1/A)(1,6) [from Eqs. (32) and (80) of [9]]. 
Comparing our solution with that given by Yamamoto [11], we note that the 
function £_(1,6, u) is missing in the integral term of Yamamoto's expression 
(20) for B{t) = £o(l,6, r). As a result, his coefficients Q(ti) and L(ji) are 
incorrect. Since £_(1,6, w) tends to 1 as 6 tends to infinity, it is clear that 
Yamamoto's formulas are a good approximation of the solution in a slab of 
large optical thickness. 



5 The case 6 — > +oo 

The relation (39) is appropriate for the computation of £o(a,b,r) when the 
values of a and 6 are such that k(a)b > 30. The dominant terms are the 
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first two terms on the right-hand side, while the third term is smallest. One 
can thus calculate l;o(a,b,T) without any (or very small) roundoff error. For 
b = +00, 0o(a, b) and tp(a, b, r) vanish, and Eqs. (39) and (25) reduce to 

£ (a, +00, r) = 1 - a (a)ip(a, +00, r), (64) 

and 

( k(a)R(a) „ r , . . . 

+f/ 1 #47'(l,r/ U )^} (65) 
2 Jo H{o>, u) u J 

respectively. Here, <So(a) — 1 — (o/2)a (a), «o(°) being the zero-order moment 
of the //-function. It is well known that this moment is (2/a)[l — \/l — a], 
so that ao(a) = y/1 — a. In Eq. (65), we have replaced 1 — ^x{o,,b,x) and 
£y(a, b, x) by their limits as b — > +00, which are 1/H(a, x) and 0, respectively. 
Using the relation (31) of [17] 

a Z 1 g(a,u) 1 i2(a) 

dw = r ~ 1 ~ ui 1 777 \ \ ' ( 66 ) 



2 7o f/"(a,u) v 7 !^ H(a,l/k(a)) 

the above expression of £0 can be written in the form 



a r 1 
"2 Jo 



exp(— t/u) \ du, (67) 



H(a, u 
which yields for a = 1 

£o(1,+oo,t) = 1. (68) 

This result means that a photon travelling in a conservative, semi-infinite 
atmosphere eventually leaves the atmosphere after multiple scattering. 



6 Numerical calculations 



The numerical evaluation of £o( a > b, r) from Eqs. (39) or (44) requires the prior 
calculation of the coefficients a (a, b) and Po(a, b), which are classical, and the 
calculation of the functions ip, tjj, ip± with the help of Eqs. (25), (28), (46), (53), 
(54). The input data are the coefficient k(a), the (classical) function H(a,x), 
and the (non-classical) functions ( + (a,b,x) and C_(a,6, x) for x G [0,1] and 
x = l/k(a). Algorithms for their computation are described in Section 6 of [9] 
and will not be repeated here. 
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In the present article, we focus on the calculation of the £o~function in the 
conservative case only, with the help of Eqs. (61) or (62). In a forthcoming 
paper, we intend to solve a problem more general than (5), with E n {r) in 
place of E 2 (t) (n > 2). Its solution £ n _ 2 (a, &, r) will be tabulated for a wide 
range of values of a, b, r and different integers n, including n = 2. Our interest 
in the conservative case is also due to its importance for stellar atmospheric 
modeling, to which we shall turn in the next section. 

We chose to calculate £o(l, b,r) from Eq. (62) instead of Eq. (61). The basic 
coefficients are (3 (l,b) and cti(l,b) + /3i(l,fe). They are tabulated in [6] and 
[7], with the first reference including useful asymptotic expressions for b — > 
and b — > +oo. More recent tables can be found in Section 9.6 of [20], which 
contains many references to previous tables. In the present work, we computed 
them accurately by means of the following formulas, given here without proof: 

^ b) = V5 [6/2 + ?NP-,o(l,6)+P-,i(l ) 6)' (69) 

(a 1 + A)(l,6) = ^P-,o(l,6), (70) 

where 

p-, n (l,b)=5 nfl -± f 1 ^^-exp(-b/u) P 4l,b,u)u n du (n>0). (71) 

2 JO -n (1, U) 

5 n fi = is the Kronecker symbol and the function p_(l,6, u) is solution to the 
following integral equation: 

ML M) = 1 - \u f expMAOMl, * ^ u (72) 

This is a Fredholm integral equation over [0, 1], with a regular kernel, which 
can be solved with great accuracy for any value of the parameters a and b [9]. 
Values of the £o-function are shown in Table 1 for a = 1 and b = 0.01, 0.1, 
0.5, 1, 10 and 100. All data are given with ten figures, which we expect to be 
significant. We have checked the 6-digit tables of Kriese and Siewert [13] and 
are able to confirm all corresponding figures given by these authors. 



7 Application to grey atmospheres in radiative equilibrium 



The function £ solves Schwarzschild's problem (5), especially for a conser- 
vative atmosphere (a = 1). This particular case was already evoked in K. 
Schwarzschild's pioneering article [1] as an interesting limiting case. Since 
Hopf's work [10], we know that the solution of an integral equation of the 
form (1) with a = 1 also models the frequency-integrated source function of a 
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Table 1 

Values of £ (1,&,t) for b = 0.01, 0.1, 0.5, 1, 10 and 100. The r-variable is less than 
6/2 since £o(l, b, b - r) = 1 - f (l, 6, r). 



6 r 


£o(l, b, t) 


b 


r 


6(1, b,r) 


0.01 


0.5126129349 


1 





0.7581464585 


0.001 


0.5097636865 




0.1 


0.6945631359 


0.002 


0.5072225665 




0.2 


0.6428723741 


0.003 


0.5047769328 




0.3 


0.5941703217 


0.004 


0.5023782905 




0.4 


0.5468089093 


0.005 


0.5000000000 




0.5 


0.5000000000 


0.1 


0.5710110310 


10 





0.9494478780 


0.01 


0.5540741012 




1 


0.8512778419 


0.02 


0.5397404273 




2 


0.7628978526 


0.03 


0.5261870381 




3 


0.6751731487 


0.04 


0.5130121656 




4 


0.5875728676 


0.05 


0.5000000000 




5 


0.5000000000 


0.5 


0.6873318619 


100 





0.9943073833 


0.05 


0.6411364778 




10 


0.8943960588 


0.1 


0.6034439404 




20 


0.7957970430 


0.15 


0.5680837612 




30 


0.6971980286 


0.2 


0.5338125198 




40 


0.5985990143 


0.25 


0.5000000000 




50 


0.5000000000 



grey atmosphere in radiative equilibrium. The conservation of energy in such 
an atmosphere reads J(b, r) = S(b, r) where J and S are the mean intensity 
and the source function, respectively, integrated over frequency. Integrating 
the formal solution of the transfer equation over the angular variable, one 
obtains [16] 

J{b,r) = J (r) + l [ b E 1 (\r-T'\)S(b,r')dT\ (73) 
I Jo 

where Jo{t) denotes the mean intensity of the direct radiation due to the ex- 
terior sources, integrated over frequency. In our problem, the exterior sources 
emit black body radiation of temperature T b across the boundary plane r = b, 
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so that J (r) = (l/2)B(T b )E 2 (b - r), where S(T 6 ) = (a/n)T b 4 (a = Stefan- 
Boltzmann constant). Accordingly, the energy condition J(b,r) = S(b,r) be- 
comes 

S(b, r) = l -B(T b )E 2 (b - T )+ l -J^E 1 {\r- r'\)S{b, r') dr', (74) 

which is an integral equation of the form (1) with a — 1. Its solution is 

S{b, r) = B(T b )Ul, b,b-r) = B(T b )[l - £ (1, b, r)], (75) 
the latter equality resulting from Eq. (42). 

From the above solution, one can derive a simple relation between the effec- 
tive temperature of the atmosphere and the temperature of the black body 
that illuminates the face r = b. The effective temperature is defined by the 
condition 

2tt ^ J o I(0,fi,u)fidfidu = aT^ f (76) 

[y = frequency variable), this quantity being the value, at r = 0, of the 
integrated radiative flux F r (r) = 2ir J +o ° J(t, fi, u)/j, dji du since there is 
no radiation incident on the face r = 0. It follows from the formal solution of 
the transfer equation that [16] 

F r {r) = 2nB(T b )E 3 (b - r) - 2n f E 2 (r - r')S(b, r') dr' 

jo 

+ 2tt jf* E 2 (r' -r)S(b,T')dr', (77) 

and condition (76) becomes 

2nB(T b )E 3 (b) + 2tt E 2 (r)S(b,r)dr = aT* ff . (78) 

Replacing in the integral S(b, r) by its second expression (75) and remarking 
that 

f E 2 {r)dr = \-E,{b), (79) 
Jo 2 

[ b E 2 (r)Co(l, b, t) dr = \[l - #,(1, 6)( ai + A)(l, 6)], (80) 
jo 2 

one obtains 

[ b E 2 (r)S(b,r)dr = B(r b )[-E 3 (b) + 6)(«i + A)(l, 6)]. (81) 

jo 2 

Replacing now nB{T b ) by aT b A , Eq. (78) finally leads to 

11/4 



T b 1 b 



T eff [/3 (1, 6)( ai + A)(l, 6)]V4 ^ _ b y 



(82) 
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The last equality in Eq. (82) follows from the relation b/3o(l, b) = {a,\— b) 
[16]. As to the Eqs. (79) and (80), only the second one is not self-evident: it 
can be found in [6], Eq. (45a). But since it is not clear to us where this equa- 
tion comes from in Ref. [6], it is proved in the Appendix. 
Eq. (82) expresses the link between the effective temperature of a grey atmo- 
sphere in radiative equilibrium and the temperature of the black body that 
irradiates it on the face r = b. It does not suppose that the atmosphere is in 
local thermodynamical equilibrium (LTE). In Fig. 1, the ratio T b /T e ff is plot- 
ted as a function of b for < b < 100. This ratio is equal to 1 for b = [since 
A)(1,0) = («i + /3i)(l,0) = 1], and it is close to (36/4) 1/4 as b — > oo [since 
«i(l,+oo) = 2/\/3 and /3i(l,+oo) = 0]. Figure 1 contains two additional 
curves (see the inlaid box), illustrating the accuracy of two expressions approx- 
imating the ratio T\,jT e ff- one is Eddington's formula Tb/T e jf = (1 + 36/4) 1 / 4 , 
which is derived e.g. in [21], p. 392, the other has been developed by one of 
us in Ref. [22] [equation preceding Eq. (39)] using a variational approach. We 
note that the accuracy of the latter approximation is better than 10~ 4 what- 
ever b > 0. 

When the atmosphere is in LTE, one can derive its temperature T(b, r) by 




20 40 60 80 

b 



Fig. 1. Tf)/T e ff versus optical thickness b. In the box: dotted and dashed curves 
represent the accuracy of two approximate expressions of T^/T e ff. Eddington's and 
Pelkowski's, respectively. 

remarking that 

S{b, r) = B[T(b, r)\ = ( ( t/tt)T 4 (6, r), (83) 
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so that, from Eq. (75) 



= Ml, b, b - r)]^ = [1 - 6(1, b, rp\ (84) 

Consequently, 

T\b,r)+T\b,b-r)=T b \ (85) 
From Eqs. (84) and (35), the surface values of the temperature are 

T fc = ' ^ ' T b = ^ °^ ' ^ ' ^ ^ 

The last equation shows that the temperature presents a discontinuity jump 
on the boundary surface t — b. This feature is well known to thermal engineers 
working on the radiative heat transfer through a gas enclosed between heated 
walls: see [6] and references cited therein. It is also known in the context of 
thermal radiation in planetary atmospheres, see for example [21] and [22]. We 
note that this temperature slip at r = b is only due to the finite character of 
the slab and has nothing to do with the hypotheses added in this section: grey 
atmosphere, radiative equilibrium, LTE. This can be seen by calculating S(b) 
from Eq. (4). Supposing B* = 0, we find S{b)/B{T b ) = (l/2)a (a,b), which is 
somewhat more general than the second Eq. (86). Our interpretation of this 
feature in a finite slab with a free surface is as follows: because of the escape 
of photons from the top boundary plane of the slab, the radiation field cannot 
be strictly isotropic within the slab, in particular very close to the boundary 
plane r = b (if b is finite). The radiation propagating in the direction of de- 
creasing r is thus anisotropic as long as it is calculated at level r < b, and it is 
"abruptly" isotropic at r = b owing to the choice of the boundary conditions. 
This originates a discontinuity in the intensity integrated over incoming direc- 
tions, hence in the mean intensity, the source function and the temperature. 
The above results are useful for the modeling of planetary and stellar atmo- 
spheres, because models are always constructed in an iterative way, and a first 
temperature distribution is necessary to start the iterations. Many modern 
codes use the Hopf distribution 

T{t) _ r 3 r _ i 1/4 



Teff ' I 

where q is Hopf 's function [23] 



{t[^ + 9W]} 1/4 , (87) 




u 



In an atmosphere code, this function is currently approximated by 2/3 what- 
ever r, and Hopf's relation becomes Eddington's law. The exact temperature 
distribution (87) was derived by Hopf [10] for a semi- infinite, grey atmosphere 
in both radiative and local thermodynamical equilibrium, with no internal 
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sources except at infinity. Equation (84) leads to its generalization to a fi- 
nite atmosphere illuminated by a black body of temperature T&. The latter 
is connected with the effective temperature of the atmosphere via Eq. (82). 
Eliminating T b from Eqs. (82) and (84), one obtains 



1-£o(1,&,t) 



T(b,r) _ f 

T eff \/5 (l,6)K + A)(l,6) 



(89) 



or, from Eq. (61), 



T(b,r) 



1 V-(lAr) 



T e// [2 (ai+A)(l, 6) 



(90) 



As b -> +oo, ^_(l,6,r) -> ^(l,+oo,r) = v^b" + g(r)] from Eqs. (65), (55) 
and (88), ai(l,6) -> «i(l) = 2/y/3, Pi(l,b) -> and Eq. (87) is retrieved. 
We thus have T(+oo,r) = T(r) and Eq. (90), with ip_(l,b, r) given by (59), 
generalizes Hopf's solution (87) to a finite atmosphere. 

We have plotted in Fig. 2 the function r -> T(b,r)/T eff for 6 = 0.01, 0.1, 1 
and +oo, and observe that T(b, r) is very close to T(r) whenever b > 1. 

1.1 




Fig. 2. T(b,r)/T e f f versus r for 6 = 0.01, 0.1, 1, +oo. 



8 Conclusion 



The problem of a slab irradiated by an isotropic radiation field is by now a 
classical one. In the present article, its solution is expressed for unspecified 
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values of a and b, without thereby excluding values of astrophysical interest 
(a ~ 1, b 1). We intend to generalize this solution to boundary conditions 
of the form = fi n (n > 0), since in the theory of stellar atmospheres 

the boundary conditions on the surface r = b are weakly anisotropic (diffusion 
approximation). 

The conservative case is particularly interesting for the modeling of planetary 
and stellar atmospheres, because it is equivalent to finding the integrated 
source function of a grey atmosphere in radiative equilibrium. It thus provides 
a good test for models which do not presuppose LTE. It is also useful when 
the LTE condition holds, since it provides the temperature distribution (90) 
that generalizes, for a finite slab, Hopf's temperature distribution (87). 



Appendix. Proof of the relation (80) 



We first prove the following lemma: for any « ^ 



a 
2 Jo 



b _ 

E 2 (r)B(a, b, r, u) dr = u[l — a (a, b)X(a, b, u) — p (a, b)Y(a, b, u)]. (Al) 



Since E 2 (t) = Jq 1 exp(— t/v) dv, the left-hand side L of the above equation is 

L —— / / exp(— t/v) dvB(a, b, r, u) dr, 
2 jo Jo 

/ / B(a, b, r, u) exp(— t/v) dr dv, 
Jo Jo 

B(a, b, 1/v, u) dv, 



a 

~2 
a 

~2Jo 

where B(a,b,l/v,u) is the finite Laplace transform, at 1/v, of the function 
r — > B(a,b,T,u). It is known (see the theorem 37.1 of [16]) that 

X(a,b,u)X(a,b,v) -Y(a,b,u)Y(a,b,v) 

B (a, b, 1/v, u = uv. A2 

v 1 ' u+v v ' 

The left-hand side of Eq. (Al) is then 

L = X(a,b,u)—u ( X(a, b, v) - — — Y(a,b, u)—u [ Y(a, b, v ) - ^ V , 
2 Jo u + v 2 Jo u + v 

= \x(a,b,u)[a (a,b) £ x (a,b,u)] 

- ^uY(a,b,u)[(3 (a,b) - -t, Y {a,b,u)} 

if we remember the definitions (13)-(15) of £x, «o and (3 . Equation (Al) 
is finally derived thanks to the classical X-equation [16] 

X(a, b, u)£x(a, b, u) — Y(a, b, u)^y(a, b, u) = X(a, b, u) — 1. (A3) 
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Integrating both members of Eq. (Al) over u from to 1 and taking into 
account Eqs. (15) and (32), one readily derives Eq. (80). 
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